arXiv: 1507.0299lv2 [cond-mat.stat-mech] 26 Oct 2015 


Rough interfaces, accurate predictions: The 
necessity of capillary modes in a minimal model of 
nanoscale hydrophobic solvation 

Suriyanarayanan Vaikuntanathan *, Grant M Rotskoff ^ Alexander Hudson^, and Phillip Geissler ^ 

* Department of Chemistry, The University of Chicago, Chicago, IL, 60637,f Biophysics Graduate Group, University of California, Berkeley, CA 94720, and f Department of 
Chemistry, University of California, Berkeley, CA 94720 

Submitted to Proceedings of the National Academy of Sciences of the United States of America 


Modern theories of the hydrophobic effect highlight its dependence 
on length scale, emphasizing the importance of interfaces in the 
vicinity of sizable hydrophobes. We recently showed that a faith¬ 
ful treatment of such nanoscale interfaces requires careful attention 
to the statistics of capillary waves, with significant quantitative im¬ 
plications for the calculation of solvation thermodynamics. Here 
we show that a coarse-grained lattice model like those of Chandler 
and coworkers, when informed by this understanding, can capture a 
broad range of hydrophobic behaviors with striking accuracy. Specif¬ 
ically, we calculate probability distributions for microscopic density 
fluctuations that agree very well with results of atomistic simula¬ 
tions, even many standard deviations from the mean, and even for 
probe volumes in highly heterogeneous environments. This accuracy 
is achieved without adjustment of free parameters, as the model is 
fully specified by well-known properties of liquid water. As examples 
of its utility, we compute the free energy profile for a solute cross¬ 
ing the air-water interface, as well as the thermodynamic cost of 
evacuating the space between extended nanoscale surfaces. These 
calculations suggest that a highly reduced model for aqueous solva¬ 
tion can enable efficient multiscale modeling of spatial organization 
driven by hydrophobic and interfacial forces. 

hydrophobic effect | lattice model | water 

Hydrophobic forces play a crucial role in biological self- 
assembly, protein folding, ion channel gating, and lipid mem¬ 
brane dynamics [1, 2, 3, 4, 5, 6 , 7, 8 , 9]. The origin and 
strength of these forces are well understood at extreme length 
scales, based on the recognition that accommodating an ideal 
volume-excluding hydrophobe in water carries the same ther¬ 
modynamic cost as evacuating solvent from the correspond¬ 
ing volume. On the scale of a small molecule like methane, 
density fluctuations that enable such evacuation are Gaussian 
distributed to a very good approximation, even far from the 
mean [10]. Linear response theories, such as Pratt-Chandler 
theory, can thus be quite accurate for assessing solvation of 
individual small hydrophobic species. 

Solvation at much larger scales is by contrast dominated 
by water’s proximity to liquid-vapor coexistence. Hydropho¬ 
bic forces involving extended substrates are shaped by the 
physics of interfaces, and are quantified by macroscopic pa¬ 
rameters like surface tension. In between these extremes, a 
rich variety of hydrophobic effects results from the combined 
importance of nearby phase coexistence and details of inter- 
molecular structure. Capturing this interplay, for instance 
near a biological macromolecule, presents a significant chal¬ 
lenge for theory. 

Lum-Chandler-Weeks (LCW) theory represents the mod¬ 
ern understanding of hydrophobic solvation, providing a con¬ 
ceptual and mathematical framework to couple interfacial 
forces with the short-wavelength density fluctuations that de¬ 
termine solvation of small molecules [10, 11]. Its physical 
perspective has inspired the development of coarse-grained 
lattice models, whose applications have revealed interesting 
and general mechanisms for the role of water in hydropho¬ 


bic self-assembly processes [2, 12, 13]. Primitive versions 
of these models, however, long appeared unable to achieve 
close quantitative agreement with atomistically detailed sim¬ 
ulations, e.g., for the probabilities of extreme number density 
fluctuations in nanometer-scale probe volumes. This short¬ 
coming motivated the construction of more elaborate mod¬ 
els, which include interactions between nonadjacent lattice 
sites ]13] and/or explicit coupling between density fluctuations 
at short and long wavelengths [12]. Quantitative accuracy im¬ 
proved as a result of these additions, but close agreement with 
detailed simulations remained elusive, despite the expanded 
set of adjustable parameters. 

The prominence of interfacial physics in this understand¬ 
ing of hydrophobic effects suggests that the quantitative suc¬ 
cess of a coarse-grained model hinges on its ability to accu¬ 
rately capture the natural shape fluctuations of a liquid-vapor 
interface ]1, 14, 15]. We recently showed that doing so with 
LCW-inspired lattice models requires closer attention to the 
statistical mechanics of capillary waves than was previously 
paid [16, 14]. These fluctuations are pronounced in molecular 
simulations, but present in lattice models only for sufficiently 
weak coupling e between lattice sites, i.e., only at tempera¬ 
tures above the roughening transition Tr [17]. In this rough 
regime the relationship between the microscopic cohesive en¬ 
ergy e and the macroscopic surface tension 7 is nontrivial. 
This previously unrecognized connection, which is essential 


Significance 

Hydrophobic effects, which play a crucial role in many chemical 
and biological processes, originate in the statistics of microscopic 
density fluctuations in liquid water. Chandler and coworkers 
have established the foundation for a simple and unified under¬ 
standing of these effects, by identifying essential aspects of wa¬ 
ter’s intermolecular structure while accounting for its proximity 
to phase coexistence. Here we show that coarse-grained mod¬ 
els based on this perspective, when constructed to respect the 
statistics of capillary waves at interfaces, can achieve remark¬ 
able agreement with results of atomistically detailed simulations. 
Highly efficient and lacking adjustable parameters, such models 
hold promise as powerful tools for studying multiscale problems 
in hydrophobic solvation. 
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Fig. 1. Density fluctuations within a cubic probe volume, of size 12A X 12A X 12A, in bulk liquid water, (a) Cross-sectional snapshot from an atomistic molecular dynamics 
simulation, showing the probe volume V in blue, (b) and (c) Probability distribution Py i^) of the number of water molecules whose center lies in V. Results are shown for 
atomistic simulations, for the LCW-inspired coarse-grained model of Eq. [5], for the conventional lattice gas of Eq. [1], and for the continuum estimate of Eq. [9]. Lattice 
models were simulated with two parameter sets consistent with the macroscopic surface tension of the air-water interface. The lower temperature case {T < parameters, 
e/T = 6, <5 = 4A) does not support capillary fluctuations. The T > Tr, parameters (e/T = 1.35, <5 = 1.84A} by contrast describe an interface whose shape fluctuations 
are consistent with capillary wave theory. Results for Eq. [5] with the T < Tr parameters, plotted in (b), are reproduced from Ref. [13]. 


for faithfully representing the spectrum of capillary waves, 
yields a lattice model parameterization that is substantially 
different than in previous work. Ref. [16] thus presented the 
first primitive LCW-inspired lattice model that fully respects 
the statistical mechanics of capillary waves. 

Here, we put the lattice model of Ref. [16] to a number of 
exacting tests, which probe its ability to accurately describe 
density fluctuations on the nanometer scales relevant to pro¬ 
tein biophysics. Despite its coarseness and lack of adjustable 
parameters, this model achieves remarkably close agreement 
with atomistic simulations, even in scenarios with strong spa¬ 
tial heterogeneity. These tests assess the importance of details 
we have omitted, such as explicit coupling between short- and 
long-wavelength density fields. 

A careful analysis of our results underscores the interplay 
of length scales accomplished by the coarse-grained model, 
highlights the importance of capillary fluctuations, and em¬ 
phasizes the special environment for solvation presented by 
extended interfaces. Our ultimate conclusion is that diverse 
hydrophobic phenomena can be captured quantitatively at a 
coarse-grained level, with minimal attention to atomic-scale 
intermolecular structure. It is sufficient to capture the cor¬ 
rect physics at extreme length scales and link them with sim¬ 
ple excluded volume constraints. We illustrate the promise 
of such models as practical tools with an application to the 
association of hydrophobic plates. 


Coarse-Graining Water 

Coarse-grained models motivated by the Lum-Chandler- 
Weeks approach separately account for density fluctuations at 
small and large length scales. As described by Chandler and 
coworkers in Refs. ]12, 13], long-wavelength variations are rep¬ 
resented on a lattice with microscopic resolution on the order 
of a molecular diameter. In the absence of solutes, walls, or 
constraints, the corresponding binary occupation variables rii, 
which indicate either vapor-like {rii = 0 ) or liquid-like (rii = 1 ) 
density in cell i, are governed by a lattice gas Hamiltonian 

H =-e'^niTij - , [1] 

(i.J> » 


where fj, is the chemical potential for cell occupation, and 
indicates summation over all pairs of nearest neigh¬ 
bor cells. We consider ambient thermodynamic conditions, at 
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which liquid is slightly more stable than vapor, /r ~ —3e + 5^, 
where 5^ is the chemical potential offset from coexistence [ 12 ] 
unless otherwise specified. As in Ref. [16], we set e/T = 1.35, 
within the narrow range that is consistent with the statisti¬ 
cal mechanics of rough interfaces yet far from criticality. The 
lattice spacing 6 — 1.84A is chosen to reproduce the experi¬ 
mentally determined surface tension 7 ™ = 72mN/nm^ of the 
air-water interface, according to the approximate relation 

= {Pef [ 2 ] 

derived in Ref. ]16]. In addition to having the correct surface 
tension, the power spectrum of interfacial height fluctuations 
of the lattice gas exhibits capillary scaling for this parameter 
set of e, 5. 

Regions that are locally liquid-like (rii = 1) additionally 
support short-wavelength density fluctuations &p(v), which 
are assumed to obey Gaussian statistics [11]. In the absence 
of constraints, these continuous fluctuations are characterized 
by the two-point correlation function 

X(r-r') = {5p(r)Ap(r')) = pt5{Y~Y')+pi(g(Y-v')-l) , [3] 

where r and r' label positions inside the liquid, pi is the macro¬ 
scopic number density of pure liquid water, and g(r) is the ra¬ 
dial distribution function [18]. We use estimates of g(r) and its 
Fourier transform obtained from experimental measurements 
by Narten and Levy [19]. 

We consider ideal hydrophobic solutes, whose influence 
on the solvent is to exclude it from a volume v. Weak, 
smoothly varying attractive interactions between solute and 
solvent amount to a small perturbation in this context. The 
effect of, e.g., dispersion forces can therefore be reasonably 
addressed using perturbation theory [12, 13]. 

We model such idealized solutes by imposing a constraint 
of solvent evacuation: The total density within v must vanish. 


E 


UiplVi + 


I dr Sp(r 

J Va 


= 0 , 


where the sum runs over all lattice cells i that intersect v, and 
Vi is the corresponding volume of intersection. 

Gaussian fluctuations in the rapidly varying field 5p(r) can 
be integrated out exactly [11], yielding an effective Hamil¬ 
tonian Hy[{ni}] for the lattice occupation variables. In the 
presence of m ideal volume-excluding solutes [12, 13], 

Hy[{ni}] = -e'^ninj-p'Y^ni + '^'^N^XTnN+ C^ , [5] 
(ij) 
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where N is an m-component column vector with elements 

. [ 6 ] 

denotes the volume of overlap between the a**' solute and 
lattice cell i, and Xin is an m x m square matrix with elements 


{Xin)a,l3 




©(r)x(r,r')e(r') . 


[7] 


Here, 0(r) = 1 if the lattice cell containing r is occupied by 
solvent and vanishes otherwise. The constant C is given by 


C = 


ln(det (27rxin)) 

max [In (det (27rxin)), Y.a 


if E.fv, > 1, 

otherwise. 


[ 8 ] 


The last two terms in Eq. [5] improve upon previous 
lattice-based models, providing a computationally tractable 
yet quantitatively accurate approximation for solute-solute in¬ 
teractions mediated by Gaussian density fluctuations in the 
surrounding solvent. Starting from a Gaussian field the¬ 
ory [11], they may be derived by applying the constraint in 
Eq. [4] separately to each solute’s excluded volume. Details 
of this derivation are included in the Supporting Information. 

Note that the coarse-grained model defined by Eq. [5] 
has essentially no free parameters—in light of Eq. [2] there 
is very little freedom in the choice of e and 5 [16]. We only 
require the surface tension of water and the bulk radial dis¬ 
tribution function g{r). The lattice model also lacks explicit 
energetic coupling between the density fields n; and Sp{r). 
Their interdependence arises strictly from the excluded vol¬ 
ume constraint expressed in Eq. [4]. 

Below we compare calculations based on this coarse¬ 
grained model with results of atomistic molecular dynamics 
simulations. Molecular simulations were performed using the 
SPC/E model of water (Methods). Several of these calcula¬ 
tions scrutinize the extreme wings of probability distributions, 
which were accessed using standard techniques of umbrella 
sampling (Methods). 


Results and discussion 

We have performed calculations that stringently assess the 
ability of an LGW-inspired lattice model to capture details of 
hydrophobic solvation at small, large, and intermediate length 
scales. We focus on characterizing and comparing the statis¬ 
tics of density fluctuations within microscopic probe volumes, 
in part because of their direct relevance to solubility. Specifi¬ 
cally, we calculate the probability Pv{N) of observing N sol¬ 
vent molecules within a probe volume of size v. Its extreme 
value determines the excess chemical potential /rex(u) of an 
ideal hydrophobe with the corresponding excluded volume, 
Atex(u) = —TlnP„(0). The behavior of Pv{N) between this 
extreme case {N = 0) and more typical values {N « pev) 
reveals much about the physical nature of fluctuations that 
might be accessed through application of external fields, so¬ 
lute attractions, or changes in thermodynamic state [2, 3]. 
With these implications in mind, we have computed Pv{N) 
over its entire meaningful range for a variety of scenarios per¬ 
tinent to solvation in complex environments. 

In addition to atomistically detailed molecular dynamics 
simulations and our LCW-inspired lattice model, we present 
results for several less sophisticated models. These simpler 
descriptions lack one or more of the physical ingredients un¬ 
derlying LOW theory, and thus shed light on their relative 
importance. Eor example, short-wavelength fluctuations can 
be straightforwardly neglected by studying the conventional 
lattice gas described by Eq. [ 1 ], which lacks biases from Gaus¬ 
sian fluctuations in 5p{r). In this case density variations 


within a probe volume can be achieved only through fluc¬ 
tuations of the binary occupation variables rn. For the pa¬ 
rameterization we have described (e/T = 1.35, 5 = 1.84A), 
which ensures T > Tr, interfaces of this lattice gas exhibit a 
spectrum of capillary modes comparable to that of a natural 
liquid-vapor interface. Fluctuations and response of these cap¬ 
illary modes may contribute significantly to solvation struc¬ 
ture and thermodynamics [4, 5, 6 [. 

To assess the importance of capillary fluctuations, we ex¬ 
amine a different parameterization of Eq. [1] (e/T = 6, 
S = 4A), for which T < Tr. This lattice gas too supports 
interfaces with the correct surface tension. Because its rough¬ 
ening transition lies above ambient temperature, however, the 
lattice gas with T < Tr lacks fluctuations in surface topog¬ 
raphy characteristic of capillary modes (i.e., the typical am¬ 
plitude of long-wavelength undulations is not proportional to 
their wavelength). Results from this model thus isolate the 
contribution of interfacial flexibility to hydrophobic effects, a 
flexibility that is also neglected in the mean field treatment of 
LOW theory [16]. 

Together, these calculations explore the interplay between 
short- and long-wavelength aspects of hydrophobicity. We find 
in general that a simple LCW-inspired lattice model can de¬ 
scribe with surprising accuracy the statistics of density fluctu¬ 
ations observed in detailed molecular simulations. This suc¬ 
cess is compromised substantially in most cases by omitting 
the effects of short-wavelength fluctuations and/or capillary 
waves, suggesting that our coarse-grained model contains a 
minimum of microscopic detail required to quantitatively cap¬ 
ture the solvation and association of nanoscale hydrophobic 
species. 


Statistics of density fluctuations in bulk water. We first exam¬ 
ine density fluctuations in the simplest aqueous environment, 
i.e., bulk liquid water. It has been well established by MD 
simulations that for small probe volumes {v < 0.5 nm^) in 
this homogeneous setting, P„(Af) has a Gaussian form well 
into its tails [10]. For larger v, low-density fluctuations are 
strongly biased by the small chemical potential difference be¬ 
tween macroscopic liquid and vapor. Pv{N) then develops an 
exponential tail, decaying much more slowly than the Gaus¬ 
sian fluctuations near {N} = ptv would suggest [13]. These 
basic features of Pv{N) are essentially built into LCW-inspired 
models, but their details can be quite sensitive to the way 
such models are constructed and parameterized. As shown 
in Fig. 1(b), the model of Eq. [5], when parameterized with 
attention to capillary fluctuations, does an excellent job repro¬ 
ducing distributions obtained from atomistic MD simulations 
for nanometer-scale cubic probe volumes, over a very wide 
range of N. 

Reversibly decreasing N in atomistic simulations from its 
average value induces formation of a small, roughly cubic 
cavity that grows to span u as Af —>■ 0. This scenario sug¬ 
gests a simple continuum estimate of Pv{N) that resolves only 
the growing interfacial area of the cavity as the probe vol¬ 
ume is evacuated. The surface area of a cubic cavity which 
accommodates an average of N water molecules in bulk is 
A = 6{N/ , which we use as an estimate of the surface 
area of the cavity that appears as water molecules evacuate 
the probe volume. Assigning the macroscopic surface tension 
7 u, as the free energy cost per unit area of the microscopic 
cavity, we obtain a prediction for the decay rate of the expo¬ 
nential tail of Pv{N), 


din P^{N) 

m 


: 4 ^ 7 ^ 


Pi 

(N) - N 


1/3 
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that agrees reasonably well with detailed simulation results. 
Lacking sensitivity to microscopic fluctuations, this estimate 
(plotted in Fig. 1(c)) unsurprisingly fails to capture the Gaus¬ 
sian character of PviN) near (N). Nor does it describe well 
the overall free energy scale associated with emptying v, erring 
by more than 70 ksT. 

Fig. 1(c) also shows results obtained from simulations of 
the lattice gas with T < Tr. Its prediction for Pv{N) closely 
resembles the simple continuum estimate of Eq. [9], accu¬ 
rately describing the low-density slope of In P„ (N) but not its 
scale or peak behavior. This similarity highlights limitations 
of lattice models at temperatures below the roughening transi¬ 
tion temperature. The deficiency of spontaneous fluctuations 
in interfacial shape begets an overly stiff response to fields or 
constraints imposed by solutes. Lattice degrees of freedom 
serve here only to coarsely determine static interfaces when 
solutes are large enough to induce drying. Lack of capillary 
modes further renders the surface tension of a lattice gas be¬ 
low the roughening transition temperature anisotropic, intro¬ 
ducing the possibility of strong lattice artifacts. In this case 
of a nanometer-scale cubic probe volume in bulk liquid, cor¬ 
respondence with the continuum estimate suggests that such 
artifacts are not substantial. 

The parameterization of Eq. [ 1 ] we have advocated, which 
does capture capillary fluctuations at the liquid-vapor inter¬ 
face, significantly improves agreement between atomistic sim¬ 
ulations and the conventional lattice gas. Plotted in Fig. 1(c), 
the T > Tr lattice gas result for Pv{N) manifests low-density 
fluctuations that are dramatically more probable than for the 
lattice gas with T < Tr. Agreement with atomistic simula¬ 
tions nonetheless remains very poor, signaling a critical role 
for short-wavelength modes even in the exponential tail of 
PviN). 

The success achieved by the full coarse-grained model of 
Eq. [5] is thus not a transparent consequence of the limiting 
behaviors motivating its form. Instead, a subtle cooperation 
of interfacial fluctuations and thermodynamics, together with 
Gaussian density statistics at the molecular scale, underlies 
its accurate prediction for Pv{N) across the entire range of 
N. 

LCW-inspired models based on the T < Tr lattice gas are 
much more difficult to reconcile with atomistic simulations. 
When parameterized with T < Tr, the unadorned form of 
Eq. [5] accurately predicts Pv{N) only near its peak, fail¬ 
ing dramatically at low N where interfacial fluctuations figure 
prominently. (See Fig. SI 1.) Ref. [13] outlines two strategies 
to address this shortcoming. Smearing out discrete interfaces 
with a numerical interpolation scheme improves predictions 
substantially, but still fails to achieve quantitative accuracy 
in the extreme tail of Pv{N). Adding as well an estimate of 
unbalanced attractive forces produces near quantitative agree¬ 
ment [13]. These elaborations, however, require introducing 
interaction potentials and adjustable parameters that are not 
clearly specified by experimental measurements [13[. Our re¬ 
sults show that greater accuracy can be achieved much more 
simply, by using lattice gas parameters that properly repre¬ 
sent the statistics of capillary fluctuations. 


Density fluctuations at a liquid-vapor interface. The inter¬ 
play of physical factors determining hydrophobic solvation 
can resolve much differently in spatially heterogeneous envi¬ 
ronments. To explore basic effects of such nonuniformity, we 
have examined microscopic density fluctuations at the inter¬ 
face between air and water. Specifically, we consider a cubic 
probe volume that straddles the plane of a macroscopic phase 
boundary, i.e., the Gibbs dividing surface between liquid and 
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vapor. The overall shape of Pv{N) in this case is similar to the 
bulk result, featuring Gaussian statistics near the peak and a 
more slowly decaying low-density tail. In this case, however, 
the thermodynamic cost of evacuation is much lower than in 
bulk, despite a similar value of {N}. 



Lattice gas, T>Tr 

0 10 20 30 40 50 60 70 

N 

Fig. 2. Density fluctuations within a cubic probe volume that straddles the in¬ 
terface between liquid water and its vapor. The probe volume, of size 12A X 12A 
X I 2 A, is centered slightly within the liquid phase, 3.67A below the Gibbs’ dividing 
surface, (a) Cross-sectional snapshot from an atomistic molecular dynamics simula¬ 
tion, showing the probe volume t; in blue, (b) Probability distribution Pv{N) of the 
number of water molecules whose center lies in V. Results are shown for atomistic 
simulations, for the LCW-inspired coarse-grained model of Eq. [5] at coexistence, 
and for the conventional lattice gas of Eq. [1] at coexistence. Lattice models were 
simulated with T > parameters (e/T = 1.35, S = 1.84A) that yield both 
the correct surface tension and capillary wave scaling. 

The LGW-inspired lattice model of Eq. [5] again matches 
atomistic simulation results very well, both near the mean of 
Pv{N) and far into its low-density wing. Our reduced descrip¬ 
tion is therefore a promising tool for assessing hydrophobic 
solvation near the liquid’s boundary. 

While the shape of Pv{N) we have determined for the 
interfacial environment resembles that of bulk liquid, the un¬ 
derlying structural fluctuations are quite different. This dif¬ 
ference is made clear by considering the simple lattice gas (in 
its higher-temperature parameterization), whose prediction is 
also plotted in Fig. 2(b). In contrast to our bulk liquid re¬ 
sults, neglecting short-wavelength density fluctuations in this 
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case effects only a modest suppression of extreme low-density 
excursions; the shape and scale of Pv{N) are in fact captured 
well by the lattice gas with T > Tr. Correspondingly, a calcu¬ 
lation based entirely on short-wavelength fluctuations, with a 
static, flat interface, fails to capture the shape of Pv{N) even 
near its peak. (See Fig. SI 1.) 

The long-wavelength density component thus dominates 
the response of the LCW-inspired model in this spatially het¬ 
erogeneous scenario, highlighting the key importance of cap¬ 
illary fluctuations at the air-water interface. Evacuation of 
a probe volume can be inexpensively achieved near a pre¬ 
existing interface by simply deforming its shape. Refs. [20] 
and [2] have also pointed to interfacial deformation as a mech¬ 
anism for extreme density fluctuations near ideal hydrophobic 
surfaces and hydrophobic biological molecules. 

The statistics of finer scale density variations have non- 
negligible quantitative impact on the predictions of Eq. [1] 
(e.g., reducing the cavitation free energy by roughly 5 ksT) 
but do not qualitatively shape the solvent response as in 
the bulk case. Underscoring the role of surface shape fluc¬ 
tuations, the lower-temperature parameterization of the lat¬ 
tice gas, which lacks capillary waves, fails profoundly to de¬ 
scribe occupation statistics for the interfacial probe volume, 
as shown in Eig. SI 1. 
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Fig. 3. Excess chemical potential of a spherical hydrophobic solute as a function 
of its perpendicular displacement z from the air-water interface, z —f oo corre¬ 
sponds to bulk vapor, z —f —oo to bulk liquid, and .2 = 0 to the Gibbs dividing 
surface between the two coexisting phases. Data are shown for atomistic simulations, 
for the LCW-inspired coarse-grained model of Eq. [5] (with T > Tr parameters 
e/T = 1.35, S = 1.84A), and for the estimate in Eq. [ 10 ] based on a completely 
quiescent interface. This solute excludes the center of each water molecule from a 
sphere of radius ft = 5A (as shown in the inset schematic). 


Association of a hydrophobic solute with the interface. The 

occupation statistics discussed above hint at thermodynamic 
driving forces that govern solvation of nonpolar species in in¬ 
terfacial environments. To make this connection explicit, we 
have computed the excess chemical potential of a spher¬ 

ical hydrophobe as a function of its perpendicular displace¬ 
ment 2 from the air-water interface (with 2 —f oo indicating 
bulk vapor and 2 —f —oo bulk liquid). Such interfacial free 
energy profiles have been determined from molecular simu¬ 
lations for a variety of solutes [4, 5, 21, 22], including small 
ions, which can exhibit a surprising tendency to adsorb to the 
liquid-vapor phase boundary [5]. While charged solutes are 
distinct in important ways from hydrophobes, the cost of cre¬ 


ating solute-sized cavities has been implicated as a key factor 
in their surface affinity [4, 5]. As a simple estimate of this 
cost, Levin et al. have assigned a fixed solvation free energy 
/disp per unit of solvent volume displaced by the solute [23]. 
For an ideally flat interface, this approximation yields a free 
energy profile, 




(flat) 

ex 


(^) 



2 (3R^ - 2^) 

TR3 


,-R<z<R, [10] 


that is antisymmetric about 2 = 0 (once fJ,ex{z) has been 
shifted by its value at 2 = 0). This estimate asserts that 
the interface influences solubility only for distances [ 2 ] smaller 
than the solute radius R. It clearly neglects the role of capil¬ 
lary fluctuations, which we have argued can be the dominant 
mode of solvent response in such scenarios. 

Atomistic simulations show that Eq. [ 10 ] poorly approx¬ 
imates the free energy profile for an ideally hydrophobic 
nanometer-scale solute. As shown in Fig. 3 for R = SA, the 
atomistic potential of mean force is highly asymmetric about 
2 = 0. On the liquid side (2 < 0), the interface’s influence ex¬ 
tends in distance 2 well beyond the solute’s radius, with /rex( 2 ) 
deviating appreciably from /rex(— 00 ) out to nearly 2 = —2R. 
This extended range strongly implicates variations in surface 
topography away from the microscopically flat geometry as¬ 
sumed by Eq. [10]. 

Our LCW-inspired lattice model, by contrast, faithfully 
captures these features of /rex(«). (See Fig. 3.) Its minimal 
ingredients thus appear sufficient to accurately describe the 
thermodynamics of accommodating volume-excluding solutes 
near the liquid’s boundary. Given the dominant role of interfa¬ 
cial softness in determining (N) for the lattice model, an ac¬ 
counting of capillary fluctuations appears essential for assess¬ 
ing the surface affinity of nonpolar solutes. At the same time, 
these results suggest that a volume-excluding hydrophobe may 
be a problematic reference system for understanding interfa¬ 
cial solvation of charged species: As a cavity near the interface 
is charged, strong response will be induced not only in solvent 
polarization, as accounted by dielectric continuum theory in 
Ref. [23], but also in interfacial shape, an aspect not addressed 
in existing theories for interfacial ion solvation. 


Density fluctuations between ideal hydrophobic plates. The 

calculations described so far demonstrate that a simple nu¬ 
merical implementation of the LCW perspective can accu¬ 
rately describe hydrophobic effects involving both microscopic 
structural response and macroscopic bistability, featuring in¬ 
terfaces that may be pre-existing or emergent. This success 
encourages use of such a model to address more complicated 
and specific situations that arise in modern biophysics and 
materials science, e.g., water flow in nanotubes [24], gating of 
transmembrane ion channels [25], or the development of ter¬ 
tiary and quaternary protein structure [26]. In each of these 
phenomena, large-scale atomistic simulations have revealed 
intriguing functional roles for hydrophobic response as solute 
configurations rearrange. Towards these frontier applications, 
we consider as a final example aqueous density fluctuations in 
a confined hydrophobic environment. 

When confined at the nanometer scale or below, water can 
exhibit physical properties markedly distinct from those of 
the bulk liquid [27]. Interactions with the containing bound¬ 
aries become a critical consideration, with hydrophilic and 
hydrophobic walls generating very different structural motifs 
and susceptibilities. Hydrophobic walls are known to greatly 
enhance density fluctuations, so that even weak external fields 
can induce nanoscale drying [3]. By poising water near a 
highly cooperative transition, such constraints can thus be 
used to engineer switching under minor perturbation, which 
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Fig. 4. Statistics of aqueous density fluctuations between nanometer-scale hydrophobic plates immersed in the liquid phase, (a) Cross-sectional snapshot from an atomistic 
molecular dynamics simulation, showing the volume-excluding hydrophobic plates in blue, (b) and (c) Probability distributions Py{N) of the number of water molecules 
whose center lies between the two plates. Results are shown for atomistic simulations and for the LCW-inspired coarse-grained model of Eq. [5] (with T > Tr, parameters 
e/T = 1.35, S = 1.84A). We also plot results for a probe volume of equivalent size in bulk water, obtained from atomistic simulations. The plates are cuboids that exclude 
the center of each water molecule from a volume with d imensions 20.24A X 5.52A X 9.2A, and are separated along their short dimension by a distance D. Plate separations 
were fixed at D = 5.52A in the simulations of panel (b), and at D = 9.2A for (c). 


in turn can sharply modulate functional behaviors like trans¬ 
port and self-assembly [1]. 

Using atomistic and coarse-grained simulations, we exam¬ 
ined a model conhnement scenario featuring two ideally hy¬ 
drophobic parallel plates (115 x 35 x 55 in size), separated 
along their short dimension by a distance D, immersed in 
liquid water, as depicted in Fig. 4(a). As in our other exam¬ 
ples, we focus on occupation statistics of a probe volume, here 
comprising the space between the two plates (a volume of size 
115 X 55 X D). Results for Pv{N) are plotted in Fig. 4 for 
two separations, D = 3S = 5.52A (b) and D = 55 = 9.2A (c). 
For comparison we also show the corresponding probability 
distributions for a probe volume of the same size and shape 
placed in homogeneous bulk water, obtained from atomistic 
simulations. 

The presence of these hydrophobic plates is not sufficient 
to induce drying for either value of D considered. The av¬ 
erage density between plates is in fact greater than in bulk 
water, due to the tendency of molecules in dense liquids to 
pack tightly against hard walls. This modest elevation of {N} 
is recapitulated by the LCW-inspired lattice model. Such an 
increase in local density is also captured by a purely Gaussian 
model that considers only short-wavelength density fluctua¬ 
tions, underscoring its origin in simple packing effects. The 
simple Gaussian model, however, considerably overestimates 
the magnitude of this shift, shown in Figs. SI 2 and SI 3— 
while these plates are not sufficiently confining to evacuate 
the probe volume with high probability, they do significantly 
enhance lattice fluctuations in their vicinity. 

The hydrophobic plates have a much stronger impact on 
the tails of Pv{N). For the smaller separation, D = 35, ex¬ 
treme low-density fluctuations are substantially more proba¬ 
ble than in the bulk case, reflecting stabilization of the dry 
state. Typical configurations do not manifest this stabiliza¬ 
tion; it would instead be apparent in the response to an ex¬ 
ternal field that disfavors occupation of the inter-plate region. 
Here, a held strength of just 0.5 T per molecule would be suffi¬ 
cient to induce drying. This profound impact of hydrophobic 
confinement on susceptibility to weak perturbations, despite 
negligible influence on typical fluctuations, has been empha¬ 
sized in previous work, and explored in detail in the context 
of protein complex formation [2, 26, 27]. 

These behaviors too are well described by the LCW- 
inspired lattice model. The slope of In Pv{N) in the low- 
density range, which largely determines the susceptibilities 
discussed above, is predicted especially accurately. The emer¬ 


gence of an inflection point, a sign of incipient bistability as 
the plates approach, is also accurately captured by the coarse¬ 
grained model. The difficulty of describing solvent-mediated 
attraction between hydrophobes using linear response theory 
has been detailed by others [28, 29]. For the nanoscale plates 
considered here, the addition of fluctuating lattice degrees of 
freedom achieves this goal with striking success. 


Conclusions 

Hydrophobic effects drive the formation of diverse assemblies 
in biological and materials systems. Our results suggest that 
the microscopic basis of these effects is thoroughly described 
by the physical perspective put forth by Lum, Ghandler, and 
Weeks. Statistics of extreme fluctuations that determine sol¬ 
vation thermodynamics can be captured with quantitative ac¬ 
curacy in a lattice model based on the LOW perspective. Do¬ 
ing so, however, requires careful attention to the softness of 
air-water interfaces, a property lacking in many previous mod¬ 
els. 

The lattice model defined by Eq. [5] appears to be truly 
minimal for this purpose. Omitting any of its contributions 
degrades the close agreement with atomistic simulations we 
have demonstrated. Moreover, its parameters are highly con¬ 
strained by basic experimental observations, namely surface 
tension, molecular pair correlations, and the spectrum of long- 
wavelength capillary waves. 

Notably absent in our model are several ingredients intro¬ 
duced in previous studies to improve agreement with detailed 
molecular simulations. We have not introduced explicit cou¬ 
pling between the rapidly and slowly varying components of 
the density field, i.e., between the lattice variables rii and 
the Gaussian held 5p(r). Others have motivated such cou¬ 
pling from the form of self-consistent equations in a mean- 
field treatment of the slowly varying density component [30], 
which manifest unbalanced attractive forces due to the rapidly 
varying component. Our omission does not, however, imply a 
lack of unbalanced attraction in the model of Eq. [5]. Direct 
interactions among the lattice variables rii in Eq. [1] are of 
course sufficient to stabilize interfaces, the primary and es¬ 
sential role of unbalanced forces in the mean field theories of 
Refs. [30] and [31]. Rather, Eq. [5] neglects the specific source 
of unbalanced attraction due to short-wavelength structure, 
a coupling whose form and strength are not transparent for 
an associated liquid like water [31]. Nor does our omission 
imply a lack of coupling between rii and 5p(r). The manip- 


6 I 


Vaikuntanathan, et al. 




ulations leading to Eq. [5] permit short-wavelength density 
fluctuations only in regions that are liquid-like (n^ = 1). By 
correlating the statistics of Sp{r) with lattice fluctuations that 
support phase coexistence, this constraint effects a potent but 
implicit interaction across length scales. Finally, in LCW- 
inspired lattice models explicit length scale coupling has the 
side effect of altering statistics of Gaussian density fluctua¬ 
tions in the vicinity of inhomogeneous lattice configurations. 
According to molecular simulations, short-wavelength fluctu¬ 
ations in liquid regions are in fact quite robust against such 
heterogeneity, particularly when liquid domains are identified 
with sensitivity to interfacial fluctuations [20, 32]. 

The lattice gas on which our model is built includes only 
nearest neighbor interactions. As emphasized in Ref. [13], the 
ground state of this model possesses unphysical degeneracies 
in the shape of closed, convex interfaces. These degeneracies 
can be removed by introducing lattice interactions between 
non-neighboring cells. In a lattice gas below its roughening 
transition temperature the impact of these additional inter¬ 
actions is dramatic, since they suppress the only affordable 
mode of interfacial shape variation. Our studies of lattice 
gases above the roughening transition temperature suggest 
that such degeneracies are much less important in the pres¬ 
ence of natural capillary fluctuations. 

Such additional couplings may be needed to further im¬ 
prove quantitative predictions, or to address more complicated 
scenarios. These goals may also require attention to molec¬ 
ular details that have not yet been incorporated into LCW- 
inspired models, for instance concerning the geometry of hy¬ 
drogen bonds, coordination statistics, or the specific form of 
interaction potentials. Pratt and coworkers have made signifi¬ 
cant advances towards understanding and quantifying the role 
of these effects in hydrophobicity [28]. Incorporating them 
into a lattice-based model poses significant challenges. 

In addition to helping establish a conceptual foundation 
for complex hydrophobic phenomena, our studies advance the 
more pragmatic goal of faithfully simulating systems that 
comprise very large numbers of water molecules. This chal¬ 
lenge limits, for example, the scale of biomolecular problems 
that can be examined by simulation without reducing the de¬ 
scription of solvent fluctuations to a gross caricature. For the 
systems we have discussed, our coarse-grained approach re¬ 
duces the computational cost of representing explicit solvent 
fluctuations by more than two orders of magnitude (relative 
to atomistic simulations), while preserving microscopic real¬ 
ism with surprising accuracy. This advantage should become 
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even more significant for very large systems, whose computa¬ 
tional burden scales exactly linearly in Eq. [5]. The model 
we have described could thus enable study of problems that 
currently lie outside the reach of atomistically detailed simu¬ 
lations. 
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Methods 

All atomistic simulations included 6,912 rigid water molecules, 
interacting through the SPC/E potential [33], in a periodi¬ 
cally replicated simulation cell with dimensions 75A x 75A 
X lOOA, and held at temperature 300 K using a Nose-Hoover 
thermostat [34, 35]. These conditions enforce phase equilib¬ 
rium, with coexisting slabs of liquid and vapor as in Ref. [20]. 
Probe volumes were positioned at least lOA away from the 
resulting liquid-vapor interface, except where indicated other¬ 
wise. Electrostatic interactions were summed using the par¬ 
ticle mesh Ewald algorithm [36]. Intramolecular constraints 
were imposed using the SHAKE algorithm [37]. Molecular dy¬ 
namics trajectories were advanced in time using the LAMMPS 
software package [38] . Probability distributions for occupation 
statistics of a probe volume were determined with umbrella 
sampling using the INDUS technique [20]. 

Coarse-grained simulations were performed with an in- 
house software package that is available upon request. These 
systems comprised 30 x 30 x 30 lattice cells, periodically repli¬ 
cated in each direction. Probability distributions of the lat¬ 
tice occupation state within a probe volume were determined 
by straightforward umbrella sampling. Statistics of the total 
density (including short-wavelength fluctuations) within such 
probe volumes were then obtained using the method outlined 
in Ref. [13]. Simulations with spherical solutes required nu¬ 
merical calculation of the overlap between cubic lattice cells 
and the solute’s excluded volume. These overlap volumes were 
computed either with Monte Carlo integration or with analyt¬ 
ical expressions detailed in the SI. Though complicated, these 
expressions significantly reduce computation time for solutes 
that move continuously in space. 
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Supporting Information: Rough interfaces, accurate predictions: The necessity of capillary modes in a minimal 
model of nanoscale hydrophobic solvation 


Integrating out Gaussian fluctuations in the presence of multiple ideal hydrophobes 

The presence of a solvent-excluding volume v introduces a constraint on the total solvent density field, 

pin(r) + Sp(r) = 0, [11] 

at each point r € v. Integrating out small-wavelength Gaussian fluctuations Sp(r) in the presence of this constraint requires 
evaluating 

Ztj = f r>(5p(r)e“'^^” 5 (pin(r)-I-(5p(r)). [12] 

The product of delta functions enforces the constraint in Eq. [11], and Hs is the unconstrained Hamiltonian for 5p{r), 

l3Hs = ^ J^J^^5p{r)x~^{r,r')5p{r'). 

A formal expression for Zv can be easily derived, but numerical evaluation is impractical due to the infinite product of 
delta functions. For this reason, prior work on coarse graining water has replaced the pointwise constraint in Eq. [12] by a 
single constraint on the average density in v, as in Eq. [MT-4] of the main text. The resulting expression for Zv is much more 
manageable and has been applied with great success to systems containing only a single solute with a simple compact geometry. 
However, for systems with many independently moving solutes that may be separated by large distances, the average constraint 
in Eq. [MT-4] introduces spurious correlations between solutes that never vanish. If v is the union of m non-intersecting regions, 

a simple solution is to constrain the average density in each region separately, i.e., to apply the constraint expressed in 
Eq. [MT-4] separately to each solute’s excluded volume. The infinite product of delta functions in Eq. [12] then becomes a 
product of m delta functions, 

5 ipin{r) + Sp{r)) 

re. . = 1 Vir6.(“) 

Using the Fourier representation of the delta function, the integral in Eq. [12] becomes 


VSp{r)e PJ I exp [ 


.(°) 


pin{r) + 6p{i 


Defining 


Nv, = 




pin{r) and <I>(r) = 


and rearranging the order of the integrals, we obtain 


n 


d'lpo' —iNnt'ilJo 


iipa, r£v^°‘\ 
0, else, 


VSp{r) exp —pHs — / 4'(r)(5p(r) . 


Evaluating the inner integral over Sp{r) results in 


Zv — Zn 


n 


dlpoi —iNa-lI^o 

-e 

27r 


exp ^ /, • 


The argument of the rightmost exponent evaluates to 


m m 




where (xin)a ,/3 is given by the double integral 

(Xin)a,/3 = / / 

With this, Zv can be expressed compactly in matrix notation, 


©(r)x(r,r')0(r'). 


Zv = Zq j ^ ^ j exp 
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where ■)/) is a column vector with elements -i/Ja, is a column vector with elements Na, and Xin is an m x m matrix with 
elements (Xin)a,/ 3 - This integral is easily evaluated to give 


Zv = Zo 


\/det(2 


'TTXi. 


exp (-iivVn'A'V 


The free energetic contribution of the small-wavelength field, relative to the unconstrained case, is then given by 


F^-Fo = -Tin ^ (in (det(2^Xin)) + fV^Xin^A^) , 


\1S] 


which are the final two terms in Eq. [MT-5] of the main text. For two solutes a and /3, the off-diagonal matrix element (xin)a,^ 
gives the coupling between Gaussian fluctuations in regions and If the solutes are separated by much more than the 
correlation length of water, this term vanishes. In the limit of infinite separation between all solutes, Xin becomes a diagonal 
matrix and the free energy in Eq. [13] reduces to a simple sum of uncoupled terms. 


T / 

Fv — Fo = — ^ ( ln(2-7r(xin)c<,a) + 
^ — 1 ^ 



Exact expression for the overlap between the lattice and spherical solute 

The Hamiltonian given in Eq. [MT-5] of the main text requires that we compute Vi, the overlap between the solute’s excluded 
volume and the cubic lattice cells. The overlap can be estimated accurately using Monte Carlo integration when the solute 
does not move over the course of the simulation. However, it is also possible to compute the overlap analytically, a pragmatic 
approach when the solute is mobile. This strategy improves accuracy and saves significant computation. 

In order to efficiently compute the overlap between cells of the lattice and spherical objects, we derived an exact formula 
using ordinary calculus. There are eight distinct cases to consider, each with several sub-cases. The calculation is a tedious, 
but straightforward exercise in enumerating the cases. There are two main, equivalent scenarios: first, if the cubic region lies 
entirely within the solute, then the overlap is simply 5^. Otherwise, the diagonal of the cubic cell, which we refer to as the line 
segment connecting {xi, yi,Zi) to (®„, y^,Zn) throughout, does not lie entirely within the sphere. Every sub-case can be rotated 
and decomposed into one of these two sub-cases. As such, knowledge of the following specific case is sufficient to compute the 
overlap in general. 

Consider a sphere of radius R centered at the origin. We will work through the case that R > S, where 5 is the coarse- 
graining length (the side length of the cubes that form the lattice). The overlap between the sphere and a cubic lattice cell 
can always be decomposed into contributions from each of the eight octants of the sphere. Without loss of generality, we will 
select the first octant x,y,z > 0, and assume that the diagonal of the cubic lattice cell lies entirely within the octant. 

With this particular case fully specified, we now carry out the integral 



dx dy dz 




dx dy dz. 


[14] 


This expression is valid when Xu,yu,Zu > R- If this is not the case, then additional integrals that account for the volume 
of the sphere in the first octant that lies outside the cubic cell must be subtracted away. These additional integrals have the 
same form as the integral in Eq. [14]. 

The integral specifying the overlap can be calculated analytically. We have not been able to simplify the expression into a 
convenient form, but the analytical result can easily be used within computer code. 
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The simplified expression for Eq. [14] is 


— ^ 2 tan ^ 


yi 

-xf- yf _ 


[ r^r^Z;- ~~ r- <"■ -(. 

- h- - 1- 

^ -i?2+a:f + zf j 


1 
2 

+ ^xi {x^ — 3R^) tan” 

zf) tan”^ 




xiyisjR? - xf - yf + xiZi^J 


I 1 f jj2 2\ , -1 f \/R^ ~ Xn — yf \ 1 /Qr)2 I 2\ , -1 ( 

+ -y,{R -y,)t.n ^ _R 2 + ^ 2 +y 2 j+^yii^R +i/0tan — 

1 /" 2 Q d 2\ i —1 \/ ~ zf \ / „9 2\ —1 / Xu\/R^- 

-Xu (*„ - 3i? ) tan I - — - 


ZR^) 

' tan ^ j 


-1 ( 

\ tan” 



Xu's/ R‘ 


yi 

-xl-yf ^ 


~ ^ \A? ^ 

/ 




[15] 


yf - XuZl\J R? - xl- zf + ^Xu\l~^\J R? -xl - zf 


^ + -1 
--tan 


/ \ 

\ H2 V-R^ _^2 -Tf y H2 -o: J ^ 

ta,n~^ r 22;;> /R2 _ -j.2 _ 2-2 _2.'r.,z,^l . 


+ i («; (3i?^ + zf)) tan ^ (2zi^R^ - 
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Fig. 5. The probability of finding N molecules in a probe volume situated at a liquid-vapor interface. Here, we show the Gaussian part of the Hamiltonian given in the main 
text, Eq. [MT-5], plotted alongside results from atomistic simulations, the coarse-grained model of Eq. [MT-5], a rough lattice gas, and a cold lattice gas at coexistence. 


12 I 


Vaikuntanathan, et al. 













a 


- 5 - 


- 10 - 


- 15 - 


- 20 - 


-25 


-30 



A -A Eq. 5, T>T^ 
Gaussian Part 
$ $ Atomistic 

Bulk Atomistic 


10 


20 


30 


40 


50 


N 


Fig. 6. The probability of finding N molecules in a probe volume between two ideal, 20.24A X 5.52A X 9.2A hydrophobic plates, separated by a distance of 5.52A. Here, 
we show the Gaussian part of the Hamiltonian given in the main text, Eq. [MT-5], plotted alongside results from atomistic simulation, the coarse-grained model of Eq. [MT-5], 
and the bulk Py[N) for a probe volume of the same size as the region between the walls. Note that the Gaussian part of the Hamiltonian overcompensates for the shift 
towards higher density. 
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Fig. 7. The probability of finding N molecules in a probe volume between two ideal, 20.24A X 5.52A X 9.2A hydrophobic plates, separated by a distance of 9.2A. Here, 
we show the Gaussian part of the Hamiltonian given in the main text, Eq. [MT-5], plotted alongside the results from atomistic simulation, the coarse-grained model of Eq. 
[MT-5], and the bulk Pv{N) for a probe volume of the same size as the region between the walls. Note that the Gaussian part of the Hamiltonian overcompensates for the 
shift towards higher density. 
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